function [sol,cal] = calibrate_abilityFixed(version, knitro_path)

%% Paths
dir_out = '../../../data/generated/matlab/Calibration/'
dir.generated = '../../../data/generated/r/';
dir.raw = '../../../data/raw/';

addpath(knitro_path)
addpath('../Densities')
addpath('../UtilityFunctions')
addpath('../../matlab')

%% Settings
% number of nodes for numerical integration
par_0.num_nodes_f = 100;
par_0.num_nodes_g = 50;
par_0.num_nodes_phi = 50;
par_0.num_nodes_bins = 100;

%% initialize various objects
% (skill distribution, utility function,participation costs)
const = tools.csv2struct([dir.raw,'constants.csv']);
tmp = tools.csv2struct([dir.generated,'mean_inc.csv']);
mean_inc = tmp.mean;
% skill distribution
par_0.dens_skill = abilityDistFixed();

% participation costs
par_0.dens_part = ParticipationNormal();

% utility function
par_util.epsilon = const.epsilon;
par_0.util = UtilityQuasiLinear(par_util);

par_0.scale_inc = 1e3; % factor for scaling income to properly
                       % calibrate tax function. (Income in the
                       % model is scaled down by this factor.)

% combined skill distribution and participation costs
dist = SkillAndParticipation(par_0);

%% Targets
% load elasticities for deciles

% decile changes from Acemoglu & Restrepo (2018), Robots and Jobs,
% in short AR (2018)
dec_eps = tools.csv2struct([dir.raw,'AcemogluRestrepo2018/dec_eps_AR2018.csv']);
targets.eps_w_K.val = dec_eps.eps;;

% frequencies and wage bins by occupation based on CPS data
freqM = tools.csv2struct([dir.generated,'freq_pareto_10_M.csv']);
freqR = tools.csv2struct([dir.generated,'freq_pareto_10_R.csv']);
freqC = tools.csv2struct([dir.generated,'freq_pareto_10_C.csv']);
targets.freqM.val = freqM.freqM;
targets.freqR.val = freqR.freqR;
targets.freqC.val = freqC.freqC;

wagebinsM = tools.csv2struct([dir.generated,'wagebins_pareto_10_M.csv']);
wagebinsR = tools.csv2struct([dir.generated,'wagebins_pareto_10_R.csv']);
wagebinsC = tools.csv2struct([dir.generated,'wagebins_pareto_10_C.csv']);
wagebins.M = wagebinsM.wagebinsM;
wagebins.R = wagebinsR.wagebinsR;
wagebins.C = wagebinsC.wagebinsC;

% modify lower and upper bounds for wage bins to have them
% consistent across occupations

w_lb = 1e-6;
w_ub = 1e4;

wagebins.M(1) = w_lb;
wagebins.R(1) = w_lb;
wagebins.C(1) = w_lb;
wagebins.all(1) = w_lb;

wagebins.M(end) = w_ub;
wagebins.R(end) = w_ub;
wagebins.C(end) = w_ub;
wagebins.all(end) = w_ub;

% participation shares by occupation and overall
part_share = tools.csv2struct([dir.generated,'part_share.csv']);
targets.part_M.val = part_share.M;
targets.part_R.val = part_share.R;
targets.part_C.val = part_share.C;
targets.participation_rate.val = part_share.all;

% prepare deciles in period 1: here need to be 10% of
% observations in each decile after the change in robots
targets.deciles_1.val = ones(8,1)*0.1;

% targets for participation change
emp_change = tools.csv2struct([dir.generated,'emp_change.csv']);
targets.part_M_change.val = emp_change.M;
targets.part_R_change.val = emp_change.R;
targets.part_C_change.val = emp_change.C;

% targets for labor income and income share
targets.avg_labor_income.val = mean_inc/par_0.scale_inc; % based on CPS
targets.transfer_lab_inc_share.val = const.transfer_share_GDP / const.labor_share; % based on

%% weight put on moments
targets.freqM.weight = 1;
targets.freqR.weight = 1;
targets.freqC.weight = 1;
targets.participation_rate.weight = 1;
targets.part_M.weight = 1;
targets.part_R.weight = 1;
targets.part_C.weight = 1;
targets.eps_w_K.weight = 5;
targets.part_M_change.weight = 0;
targets.part_R_change.weight = 1e4;
targets.part_C_change.weight = 0;
targets.avg_labor_income.weight = 1;
targets.transfer_lab_inc_share.weight = 1;


%% initialize calibration object
varnames = {'Y_M_0','Y_R_0','Y_C_0',...
    'Y_M_1','Y_R_1','Y_C_1',...
    'muM','muR','muC',...
    'sigmaM','sigmaR','sigmaC',...
    'aM','aR','aC',...
    'transfer','shareM','shareR','shareC',...
    'mu_part_M',...
    'sigma_part_M','xi','tau_hsv','lambda_hsv'};


cal = ...
    CalibrateKrusellWithRobotsEconomy(dist, ...
    targets, wagebins, ...
    varnames,w_lb,w_ub, ...
    'num_nodes_f', ...
    par_0.num_nodes_f, ...
    'num_nodes_g',par_0.num_nodes_g, ...
    'num_nodes_bins',par_0.num_nodes_bins);

%% initial point and bounds
smallNum = 1e-4;
bigNum = 1e6;

% initial point
inp.Y_M_0 = 8.162223691;
inp.Y_R_0 = 9.184529177;
inp.Y_C_0 = 14.39994356;
inp.Y_M_1 = 10.55253833;
inp.Y_R_1 = 8.471985238;
inp.Y_C_1 = 14.95999003;
inp.muM = 0;
inp.muR = 0;
inp.muC = 0;
inp.sigmaM = 0.331176365;
inp.sigmaR = 0.327738986;
inp.sigmaC = 0.369296071;
inp.aM = 4.854743237;
inp.aR = 2.371641013;
inp.aC = 2.432201324;
inp.transfer = 4.807022257;
inp.shareM = 0.132;
inp.shareR = 0.574;
inp.shareC = 0.294;
inp.mu_part_M = -17.80318728;
inp.sigma_part_M = 40.82353757;
inp.xi = 0.578473012;
inp.tau_hsv = const.tau_hsv;
inp.lambda_hsv = const.lambda_hsv;

% lower bounds
lb.Y_M_0 = 0.1;
lb.Y_R_0 = 0.1;
lb.Y_C_0 = 0.1;
lb.Y_M_1 = 0.1;
lb.Y_R_1 = 0.1;
lb.Y_C_1 = 0.1;

lb.muM = 0;
lb.muR = 0;
lb.muC = 0;

lb.sigmaM = 0.01;
lb.sigmaR = 0.01;
lb.sigmaC = 0.01;

lb.aM = 0.1;
lb.aR = 0.1;
lb.aC = 0.1;

lb.shareM = inp.shareM;
lb.shareR = inp.shareR;
lb.shareC = inp.shareC;

lb.transfer = 0;

lb.mu_part_M = -bigNum;
lb.sigma_part_M = smallNum;
lb.xi = 1e-6;

lb.tau_hsv = inp.tau_hsv;
lb.lambda_hsv = inp.lambda_hsv;

% upper bounds
ub.Y_M_0 = 50;
ub.Y_R_0 = 50;
ub.Y_C_0 = 50;
ub.Y_M_1 = 50;
ub.Y_R_1 = 50;
ub.Y_C_1 = 50;

ub.muM = 0;
ub.muR = 0;
ub.muC = 0;

ub.sigmaM = 10;
ub.sigmaR = 10;
ub.sigmaC = 10;

ub.aM = 20;
ub.aR = 20;
ub.aC = 20;

ub.shareM = inp.shareM;
ub.shareR = inp.shareR;
ub.shareC = inp.shareC;

ub.transfer = 1e6;

ub.mu_part_M = bigNum;
ub.sigma_part_M = bigNum;
ub.xi = 1e10;
ub.tau_hsv = inp.tau_hsv;
ub.lambda_hsv = inp.lambda_hsv;

%% Calibrate
cal.moments_skill_prices(inp)
sol = cal.calibrate_skill_dist_cons(inp,lb,ub) %actual calibration (cons stands for constrained)
mom = cal.moments_skill_prices(sol)

%% Save to file
% means are set to zero
sol.muM = 0;
sol.muR = 0;
sol.muC = 0;

% previously had occupation-specific participation costs, now no
% longer, hence just assign same values mu_part_M and sigma_part_M
% to all
sol.mu_part_R = sol.mu_part_M;
sol.mu_part_C = sol.mu_part_M;

sol.sigma_part_R = sol.sigma_part_M;
sol.sigma_part_C = sol.sigma_part_M;

% also save scale_inc
sol.scale_inc = par_0.scale_inc;

save([dir_out,version,'_results.mat'],'cal','sol')

tools.struct2csv(sol,[dir_out, version, '_abilityDistThreeFixC_sol.csv']);

fields = {'freqM','freqR','freqC','deciles_0','deciles_1', ...
    'eps_w_K'};
mom = rmfield(mom,fields);
tools.struct2csv(mom,[dir_out, version, '_abilityDistThreeFixC_mom.csv']);
end